UNLIMITED 


RSRE 

MEMORANDUM  No.  4409 


gj  ROYAL  SIGNALS  &  RADAR 
f  ESTABLISHMENT 

: .  Q 

;  < 


THE  QR-DECOMPOSITION  BASED  LEAST-SQUARES 
LATTICE  ALGORITHM  FOR  ADAPTIVE  FILTERING 


Authors:  I K  Proudler,  J  Q  McWhlrter 


an 

o 

KT 


6 


UJ 

he 
m 
t c 


PROCUREMENT  EXECUTIVE, 
MINISTRY  OF  DEFENCE, 
RSRE  MALVERN, 
WORCS. 


sinTfirOTirin 


i  iminm'-Ti 

"t;  ,mmt  —  .. 

for  pibfle 

UUtJbt&Oa  ttii&iiSW 


UNLIMITED 


CONDITIONS  OF  RELEASE 

0083098  BR  115234 


. . . . .  DRICU 

COPYRIGHT  (c) 

1988 

CONTROLLER 
HMSO  LONDON 


. . .  DRICY 

Reports  quoted  are  not  necessarily  available  to  members  of  the  public  or  to  commercial 
organisations. 


ROYAL  SIGNALS  AND  RADAR  ESTABLISHMENT 


Memorandum  4409 

TITLE:  THE  QR  DECOMPOSITION  BASED  LEAST-SQUARES  LATTICE 

ALGORITHM  FOR  ADAPTIVE  FILTERING 

AUTHORS :  I  K  Proudler  and  J  G  McWhirter. 

DATE:  July  1990 


SUMMARY 

We  derive,  from  first  principles,  the  least  squares  lattice  algorithm  for  adaptive  filtering  based  on 
the  QR  decomposition  (QRD).  In  common  with  other  lattice  algorithms  for  adaptive  filtering,  this  algo¬ 
rithm  only  requires  0(p )  operations  for  the  solution  of  a  p-th  order  problem.  The  algorithm  has  as  its 
root  the  QRD-based  recursive  least  squares  minimisation  algorithm  and  hence  is  expected  to  have  su¬ 
perior  numerical  properties  when  compared  with  other  fast  algorithms. 

This  algorithm  contains  within  it  the  QRD-based  algorithm  for  solving  the  least  squares  linear  pre¬ 
diction  problem.  These  algorithms  are  presented  in  two  forms:  one  that  involves  talcing  square-roots  and 
one  that  does  not.  Some  computer  simulations  of  a  channel  equaliser,  using  finite-precision  arithmetic, 
are  presented  in  which  the  lattice  algorithms  are  compared  to  the  more  established  triangular  systolic 
array  ones. 

The  relationship  between  the  QRD-based  lattice  algorithm  and  other  least  squares  lanice  algo¬ 
rithms  is  briefly  discussed.  Various  extensions  to  this  work  are  discussed  including  the  multi-channel 
QRD-based  adaptive  filtering  algorithm  that  can  be  used  for  wide-band  beamforming. 
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1.  INTRODUCTION 


Least  squares  minimisation  can  be  applied  to  a  wide  range  of  digital  signal  processing  problems 
including  that  of  adaptive  filtering.  The  adaptive  filtering  problem  is,  however,  special  in  the  sense  that 
the  “data  matrix"  (Xp(n)  of  equation  (5))  has  a  Toeplitz  structure  so  that  each  row  of  the  data  matrix 
contains  only  one  new  datum  when  compared  to  the  previous  row.  Various  algorithms  [8]  have  been 
devised  that  take  advantage  of  this  redundancy  to  reduce  the  computational  load,  for  a  p-th  order  filter, 
from  Ofp2)  to  O(p).  Unfortunately  these  fast  algorithms  are  not  explicitly  well-conditioned  and  some 
tend  to  be  numerically  unstable. 

It  is  possible,  however,  to  solve  a  least  squares  minimisation  problem  using  the  technique  of  QR 
decomposition[8].  This  technique  has  the  advantage  that  it  operates  on  the  data  matrix  directly,  rather 
than  the  corresponding  covariance  matrix,  and  involves  only  orthogonal  rotations,  which  are  numeri¬ 
cally  well-conditioned.  The  recursive  QRD  algorithm  can  be  implemented  on  the  triangular  systolic  ar¬ 
ray  as  devised  by  Gentleman  and  Kung[5]  and  subsequently  modified  by  McWhirter[14].  This  algo¬ 
rithm  solves  a  general  recursive  least  squares  minimisation  problem  and  will  require  0(p2)  operations 
to  generate  the  solution  to  a  p-th  order  adaptive  filter  problem.  Extensive  computer  simulations  of  this 
more  general  algorithm[28]  have  shown  the  QRD-based  least  squares  minimisation  algorithm  to  have 
excellent  numerical  properties.  The  possibility  of  producing  an  O(p)  QRD-based  algorithm  for  the  spe¬ 
cial  case  of  adaptive  filtering  has  thus  long  been  of  interest. 

Here  we  present  a  full  derivation  of  the  recently  derived  QRD-based  least  squares  lattice  algo¬ 
rithm!  16]  starting  from  the  Ofp2)  QRD-based  algorithm  and  incorporating  directly  the  time-shift  redun¬ 
dancy  found  in  an  adaptive  filtering  problem.  The  derivation  presented  here  owes  much  to  the  work  of 
Cioffi[3).  Recently  he  showed  how  to  lake  advantage  of  the  time-shift  redundancy  that  is  present  in  the 
problem  of  the  linear  prediction  of  time  series  data  and  produced  a  QRD-based  “fast  Kalman"  algorithm 
(see  [  1]  or  [  16]  for  alternative,  simpler  derivations).  This  algorithm  can  recursively  update,  from  one 
time  instant  to  the  next,  the  solution  to  a  p-th  order  linear  prediction  problem  in  O(p)  operations.  How¬ 
ever,  unlike  other  fast  Kalman  algorithms,  the  QRD-based  fast  Kalman  algorithm  also  produces  the  so¬ 
lution  to  all  lower  order  problems  as  well  (see  [19]). 

The  QRD-based  lattice  algorithm,  like  all  lattice  algorithms,  is  designed  to  solve  the  linear  predic¬ 
tion  problem  recursively  in  time  and  order,  again  in  Ofp)  operations:  thus  it.  too,  solves  all  lower  order 
problems.  The  two  classes  of  fast  QRD-based  algorithms  are  different  however.  The  fast  Kalman  algo¬ 
rithms  are  completely  different  in  structure  from  the  QRD-based  lattice  algorithms  (see  Slock  [24]  for 
more  discussion).  In  particular  the  fast  Kalman  algorithms  have  a  smaller  operation  count  than  the  lat¬ 
tice  algorithms  (although  both  are  linear  in  the  problem  order).  The  level  at  which  the  two  different  al¬ 
gorithms  can  be  pipelined  is  also  different  -  the  lattice  algorithms  having  a  higher  degree  of  concurren¬ 
cy.  It  is  also  worth  noting  that  the  fast  Kalman  algorithms  implicitly  have  a  downdating  step  which  gives 
cause  for  concern,  from  a  numerical  point  of  view,  although  no  problems  have  been  observed  in  any 
simulations  done  to  date. 

In  recent  months,  the  Ofp)  QRD-based  least  squares  lattice  algorithm  has  also  been  derived,  inde¬ 
pendently,  by  two  other  groups:  Ling[13]  and  Yang  and  Bdhme[29],  Both  of  these  derivations  begin 
from  a  previously  known  (non-onhogonal1)  algorithm  and  by  a  series  of  transformations  arrive  at  one 
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with  only  orthogonal  operations.  As  such  these  derivations  provide  an  interesting  insight  into  the  prob¬ 
lem  and  are  complementary  to  that  presented  here. 

Vang  and  BOhme  bring  together  the  woric  of  Lewis[9]  and  McWhirter[14],  Lewis  began  with  the 
standard  (covariance  domain)  multi-channel  least  squares  lattice  equations  and,  driven  firstly  by  the  de¬ 
sire  not  to  explicitly  invert  the  covariance  matrices  and  secondly  to  perform  all  matrix  computations  in 
a  ‘numerically  sound”  way,  reformulated  part  of  the  algorithm  (the  calculation  of  the  reflection  coeffi¬ 
cients)  by  the  use  of  the  matrix  inversion  lemma  and  QR  decomposition.  As  the  bulk  of  the  calculation 
is  exactly  the  computation  of  the  reflection  coefficients,  Lewis  proceeded  no  further  with  this  re-formu- 
lation  and  apparently  failed  to  notice  that  the  “non-onhogonal”  part  of  his  algorithm  is  in  fact  redundant 
Following  the  work  of  McWhirter,  Yang  and  BOhme  noticed  that  the  adaptive  filtering  residuals  can  be 
found  directly:  this  observation  results  in  the  construction  of  a  purely  orthogonal  algorithm  (see  Figure 
1). 
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Figure  1  Relationship  to  Covariance  Domain  Lattice  Algorithm  and  RMGS 


The  other  derivation,  by  Ling,  relies  on  a  well  known  equivalence!  10]  between  two  general  least 
squares  minimisation  algorithms:  the  0(p1 2)  QRD-based  algorithm  and  the  recursive  modified  Gram- 
Schmidt  (RMGS)  algorithm!  1 1].  Strictly  speaking,  the  equivalence  is  between  a  modification  of  the 
RMGS  algorithm  and  a  form  of  the  QRD-based  algorithm  that  does  not  require  the  square-root  opera¬ 
tion2.  Ling  and  Proakisf  12]  have  shown  how  to  create  an  0(p)  lattice  algorithm  for  solving  the  adaptive 
filtering  problem  by  taking  advantage  of  the  time-shift  redundancy  within  the  RMGS  approach:  the  so 
called  “error  feedback”  lattice  algorithm.  This  is  arguably  the  most  well-behaved  of  the  least  squares 
lattice  algorithms  presented  to  date.  Using  the  above  equivalence.  Ling  [13]  has  recently  shown  that  the 
error  feedback  lattice  algorithm  can  be  reinterpreted  as  being  a  square-toct-free  form  of  an  associated 
algorithm  that  consists  only  of  orthogonal  rotations.  Comparison  with  the  algorithm  derived  here  shows 
that  Ling's  associated  orthogonal  algorithm  and  that  of  Yang  and  BOhme,  when  specialised  to  the  single 

1.  We  use  the  phrase  “orthogonal  algorithm"  to  mean  one  that  generates  the  required  solution  exclusively 
by  the  application  of  orthogonal  transformations,  such  as  Givens  rotations,  to  the  input  data. 

2.  See  section  5  for  details  of  the  square-root-free  version  of  the  QRD-based  algorithm. 


2 


channel  case,  are  indeed  identical  and  are  equivalent  to  the  QRD-based  least  squares  lattice  algorithm. 

Another  approach  to  fast  orthogonal  adaptive  filtering  algorithms  has  been  presented  by  Regalia 
and  Bellanger  [20],  Based  on  the  work  of  Cioffi[3],  Regalia  and  Bellanger  realised  that  certain  quanti¬ 
ties  in  the  QRD-based  fast  Kalman  algorithm  were  the  same  as  those  in  a  conventional  (covariance  do¬ 
main)  lattice  algorithm.  In  particular,  the  identification  of  the  reflection  coefficients  as  the  sines  of  cer¬ 
tain  rotation  angles  led  them  to  develop  an  alternative,  Kalman-type  algorithm  for  solving  the  linear  pre¬ 
diction  problem.  Regalia  has  shown  theoretically  [21)  that  his  structure  is  stable.  However  it  is  not  as 
yet  clear  whether  the  same  analysis  will  work  for  the  other  fast  QRD-based  adaptive  filter  algorithms. 

As  well  as  presenting  the  QRD-based  least-squares  lattice  algorithm  for  linear  prediction,  the  ex¬ 
tension  of  this  technique  to  the  solution  of  the  adaptive  filtering  problem  is  also  given  in  this  memoran¬ 
dum.  The  resulting  algorithm  has  a  lattice-ladder  structure.  In  common  with  Cioffi's  original  formula¬ 
tion  of  the  QRD-based  fast  Kalman  algorithm,  the  lattice-ladder  algorithm  presented  here  operates  on 
pre-windowed  data  (i.e.  all  data  before  the  first  time  instant  is  assumed  to  be  zero).  The  extension  of  this 
work  to  the  multi-channel  case  (wide-band  beamforming)  is  relatively  straightforward  and  is  briefly  dis¬ 
cussed  in  section  6. 

We  begin,  in  section  3,  by  reviewing  the  mathematics  of  the  solution  to  an  adaptive  filtering  prob¬ 
lem  using  the  method  of  QR  decomposition.  The  key  to  the  fast  adaptive  filtering  algorithm  is  the  de¬ 
velopment  of  a  fast  (lattice)  algorithm  for  the  solution  of  an  associated  linear  prediction  problem.  In  sec¬ 
tion  4  the  connection  between  these  two,  related,  problems  is  outlined  and  the  O(p)  solution  to  the  linear 
prediction  problem,  and  hence  the  adaptive  filtering  problem,  is  developed.  Section  5  discusses  various 
aspects  involved  in  the  implementation  of  the  lattice  algorithm.  The  derivation  of  the  multi-channel  lat¬ 
tice  algorithm  is  sketched  out  in  section  6  and  the  results  of  some  computer  simulations  are  presented 
in  section  7.  The  algorithm  is  given  explicitly,  in  the  appendix,  in  two  forms:  one  that  involves  the 
square-root  operation  and  one  that  does  not. 
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2.  NOTATION 


11*11  L2  norm. 

Q  a  vector  of  zeros3. 

O  a  matrix  of  zeros. 

an  n-dimensional  vector  full  of  zeros  except  for  the  n-th  component  which  is 
unity  (pinning  vector). 

Q  an  orthogonal  rotation  matrix. 


Q 


time  recursive  increment  of  Q.  E.g.  Q(n)  =  Q(n) 


Q(n-  1)  £> 
£>l  1 


an  upper  triangular  matrix, 
a  data  matrix, 
a  reference  vector. 

the  two  components  of  the  rotated  reference  vector:  Qy. 
the  two  components  of  the  right-hand  column  of  the  matrix  Qp.  Alternatively, 
the  two  components  of  the  rotated  pinning  vector:  Qn 
a  least  squares  residual, 
a  routed,  or  angle-normalised  residual, 
the  exponential  weighting  factor, 
the  lower  right-hand  element  of  the  matrix  Qp. 
the  sum  of  the  squared  prediction  errors. 

Subscript  T'  indicates  a  quantity  connected  with  a  forward  linear  prediction  problem  or  its 
solution. 

Subscript  "b"  indicates  a  quantity  connected  with  a  backward  linear  prediction  problem  or  its 
solution. 

Subscript  "p”  indicates  a  quantity  connected  with  a  p-th  order  problem  or  its  solution  (simi¬ 
larly  subscripts  "p-1",  "p+1"  etc.). 

Hence,  for  example: 


R 

X 

i 

U,  y 

&P-&P 

e 

a 

P 

TP 

e 


p  j(n-l)  is  the  top  component  of  the  routed  reference  vector  y^  p  j(n- 1 )  associated  with 

the  (p-l)st  order  backward  linear  prediction  problem  at  time  (n-1). 

Several  roution  matrices  are  introduced  in  this  memorandum  that,  clearly,  have  to  be  distinguished 
from  each  other.  There  is,  however,  no  obvious  choice  as  to  how  they  should  be  labelled.  Following  the 
convention  used  for  labelling  of  the  reflection  coefficients  in  the  covariance-domain  least  squares  lattice 
algorithm,  we  choose  to  label  these  roution  matrices  according  to  the  problem  to  which  they  relate.  For 
example  the  matrix  Qf  p(n)  is  a  roution  matrix  used  in  the  calculation  of  the  solution  to  the  p-th  order 
forward  linear  prediction  problem  at  time  n:  despite  the  fact  that  it  operates  on  the  vector 
Yb.a-,(n-l)! _ 

3.  The  dimensionality  of  the  all  zero  vectors  and  matrices  used  in  this  memorandum  should  be  obvious 
from  the  context  in  which  they  are  used. 
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3.  ADAPTIVE  FILTERING 


In  order  to  simplify  the  analysis  we  consider  only  real  signals.  The  extension  to  the  complex  case 
is  straight  forward  and  indeed  the  algorithms  presented  in  the  appendix  are  for  complex  signals. 

In  a  least  squares  adaptive  filtering  problem,  a  set  of  p  (say)  weights,  ffip(n)  = 

(o)p0,(n),  £Op^(n), COpP-1*(n)] ',  is  to  be  found,  at  time  n,  that  minimises  that  the  sum  of  the  dif¬ 
ferences,  squared,  between  a  reference  signal  y(n)  and  a  linear  combination  of  p  samples  from  a  data 
time  series  x(n-i)  (0  £  i  £  p-1).  This  decision  is  to  be  based  on  all  the  data  accumulated  so  far. 
Specifically,  the  measure  Ep(n)  is  to  be  minimised,  where: 

Ep(n)=IIB(n)ep(n)ll  (1) 


p"-1  0  . 

..  0 

B(n)  = 

0  pn~2. 

..  0 

(2) 

0  0  . 

..  1 

0<p<l 

(3) 

£p(n)  = 

Xp(n)cop(n)  +  y(n) 

(4) 

Xp(n)  = 


x(I) 

x(2) 


x(0) 

x(l) 


x(2-p) 
x(3  -  p) 


x(n)  x(n  -  1 )  ...  x(n  -  p  +  1  )J 


(5) 


i(n)=[y(l) . y(n)]1  (6) 

The  diagonal  matrix  B(n)  represents  an  exponential  "forgetting"  factor  that  allows  the  algorithm 
to  work  with  quasi-stationary  signals.  Note  that  there  are  effectively  two  time  indices  present  in  equa¬ 
tion  (4)  because  it  is  necessary  to  distinguish  between  the  components  of  the  error  vector  gp(n)  and  the 
amount  of  data  used  to  calculate  the  coefficients  (and  hence  the  time  index  forap(n)l.  Specifically  we 
have 


fp(n)  =  (ep(l,n),ep(2,n),  ...,ep(n,  n)] 1 


(7) 


p-  1 

where  ep(m,  n)  =  y(m)+  £  o)(,)(n)x(m  -  i) 

i  =  0 


(8) 
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is  the  error  in  estimating  the  datum  y(m)  as  a  linear  combination  of  x(m-i)  0Ki<p-l  when  using  the  op¬ 
timum  coefficients  calculated  at  time  n.  There  are  two  quantities  of  panicular  interest:  the  a-posteriori 
error  (equation  (9))  which  uses  the  most  up  to  date  estimate  of  y(n)  and  the  a-priori  error  (equation  (10)) 
which  is  the  estimation  error  for  y(n)  based  on  the  coefficients  calculated  at  the  previous  time. 

p- 1 

ep(n,  n)  =  y(n)+  £  ©^(njxfo-i)  (9) 

i  =  0 


P-  1 


ep(n,  n  -  1 )  =  y(n)+  £  cop  (n-  l)x(n-i) 
i  =  o 


(10) 


The  solution  of  least  squares  minimisation  problems  via  QR  decomposition  is  now  well  established 
[8]  and  requires  the  determination  of  an  orthogonal  matrix  Qp(n)  that  transforms  the  matrix  B(n)  >-p<n) 
into  an  upper  triangular  form  (as  indicated  by  the  shading).  Let 


Qp(n)  B(n)  Xp(n)  = 


'Rp(n) 

O 


so  that,  from  equations  (4)  and  (11), 

Up(n)‘ 

_Vn)_ 


Qp(n)B(n)ep(n)  = 


Rp(n) 


w(n)  + 


(11) 


(12) 


where 


yp(n) 

_Vn). 


E  Qp(n)  B(n)  y(n) 


(13) 


Note  that  Bpfn)  is  defined  to  be  a  p-dimensional  vector  so  that  the  partitioning  in  equation  (12)  is  con¬ 
sistent.  Due  to  the  orthogonal  nature  of  Qp(n)  it  is  clear  that 


Ep(n)  =  IIB(n)  Sp(n)ll  =  IIQp(n)  B(n)  gp(n)ll  = 


Rp(n) 

O 


fc>(n) 


np(n)l 

,yp(ny 


(14) 


and  by  inspection  it  can  be  seen  that  the  least  squares  solution  is  obtained  when 

Rp(n)ffl(n)  +  Upfn)  =  2 


(15) 


and  that 


n$j{Ep(n)}  =  ||yp(n)|| 


(16) 
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In  principle  equation  (15)  can  be  solved  by  back  substitution,  since  Rp(n)  is  a  triangular  matrix, 
and  the  optimum  coefficient  vector ia(n)  found.  It  is  often  the  case,  and  adaptive  filtering  is  a  prime  ex¬ 
ample,  that  the  a-posteriori  residual  (equation  (9))  is  explicitly  required4.  The  reason  for  this  is  the  fact 
that  the  a-posteriori  residual  represents  that  part  of  the  reference  signal  uncorrelated  with  the  data  sig¬ 
nals  x(n-i)  (OSiSp-1)  and  as  such  is  the  "filtered”  signal.  Dearly  the  a-posteriori  residual  can  be  calcu¬ 
lated,  oncecfip(n)  is  known,  as  the  last  component  of  £p(n)  (equation  (7)).  It  has  been  shown  [14],  how¬ 
ever,  that  the  adaptive  filtering  residual  can  be  obtained  directly  from  the  product  of  two  quantities  nat¬ 
urally  present  in  the  QRD-based  algorithm: 


ep(n.n)  =  Yp(n)ap(n) 

(17) 

where 

Yp(n)=2t'  Qp(n)  2n 

(18) 

Vn)=  tPVn) 

(19) 

and 

5n=[° . 0.1]' 

(20) 

is  an  n-dimcnsional  vector.  Vectors  of  the  form  shown  on  equation  (20)  consisting  of  zeros  except  for 
the  final  element,  which  is  unity,  play  an  important  rfile  in  the  following  mathematics.  These  vectors  arc 
often  referred  to  as  "pinning"  vectors.  In  equations  (18)  and  (19),  the  pmning  vectors  merely  serve  to 
select  the  lower  right-hand  element  of  Qp(n)  and  the  last  element  of  \^(n>  respectively. 

The  ability  to  calculate  the  adaptive  filtering  residual  directly,  using  equation  (17),  is  an  important 
result  since  although  inverting  the  matrix  Rp(n)  is  relatively  easy  (it  is  triangular  and  therefore  can  be 
inverted  by  back-substitution)  this  process  is  potentially  numerically  unstable  whereas  the  two  quanti¬ 
ties  yp(n)  and  ap(n)  are  always  well  defined.  It  can  be  shown’  ’]  that  yp(n)  is  the  square-root  of  the 
maximum  likelihood  factor  of  more  conventional  algorithms,  .c  follow  Ling  [13]  in  referring  toap(n) 
as  the  "angle  normalised”  residual. 

It  is  worth  noting,  at  this  point,  that  there  arc  two  methods  for  calculating  the  orthogonal  matrix  0 
via  Givens  (planar)  rotations  and  via  the  Householder  transformation  (a  reflection).  The  Householder 
technique  leads  to  a  lower  operation  count;  however,  unlike  the  Givens  approach  it  docs  not  admit  a 
time  recursive  implementation,  which  is  often  preferred5. 

It  can  not  be  stressed  too  much  that  once  Qp(n)  has  been  found  the  problem  has  effectively  been 

4.  Note  that  it  is  also  possible  to  extract  the  coefficient  vector  from  the  QRD-based  algorithm  in  a  systolic 
fashion:  see  section  10.4.1. 

5.  Recently,  Cioffi  (4]  and  Slock  [23]  have  published  recursive  QRD-based  "fast  Kalman"  algorithms 
that  involve  Householder  transformations,  but  the  steps  that  involve  time  updating  are  operations  on 
2-dimensional  vectors  and  as  such  the  required  Householder  transformations  are  equivalent  to  Givens 
rotations.  There  does,  however,  appear  to  be  some  controversy  as  to  the  validity  of  these  algorithms 
see  [23], 
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solved.  Knowledge  of  Qp(n)  means  that  yp(n)  is  known  and  also  allows  ap(n)  to  be  calculated  and  thus 

the  least  squares  residual  may  then  be  found.  The  development  of  the  fast  QRD  algorithm  for  adaptive 
filtering  is  based,  almost  entirely,  on  the  principle  of  constructing  partially  triangularised  matrices  from 
known  quantities  and  then  finding  a  set  of  rotations  to  complete  the  process.  This  principle  is  also  a  key 
element  in  the  derivation  of  the  well  known  [5]  time-recursive  version  of  the  algorithm  outlined  above. 
Here  the  solution  at  time  n,  along  with  the  new  input  data  for  time  (n+1 ),  is  used  to  simplify  the  calcu¬ 
lation  of  the  solution  at  time  (n+1).  The  approach  may  be  summarised  as  follows:  since 


Xp(n+1)  = 


Xp(n) 

x(n  +  l)...x(n  -  p  +  2) 


(21) 


it  follows,  from  equation  (11),  that 


Qp(n)  9 

pRp(n) 

.  9l  1 

B(n+1)  Xp(rn-l)  = 

0 

x(n+  l)...x(n-p  +  2) 

(22) 


Note  that  the  matrix  on  the  right-hand  side  in  equation  (22)  is  very  nearly  triangular  and  composed 
entirely  of  known  quantities.  Thus  the  actual  application  of  the  rotations  specified  in  equation  (22)  can 
be  circumvented  by  the  direct  construction  of  the  partially  triangularised  matrix  shown  in  equation  (22). 
To  complete  the  triangularisation,  define  Qp(n+1)  to  be  that  set  of  rotations  that  annihilates  the  new 
data  samples  by  rotating  them  into  the  matrix  PRp(n).  In  which  case 


Qp(n+1) 


Qp(n)  p 

Rp(o  +  l) 

-  £>'  L 

B(n+1)  Xp(n+1)  = 

L  0 

(23) 


As  only  the  rotations  in  Qp(n+1)  have  to  be  constructed  (as  a  series  of  Givens  rotations[5)),  the 
computational  load  is  reduced  from  O(np)  to  0(p2).  Note  that,  from  equations  (22)  and  (23), 


Qp(n+l)=Op(n+l) 


Qp(ri)  P 
.  £>l  11 


(24) 


so  that 


=  ®n+l 


(25) 
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Dp(n  +  1 ) 

=  Qp(n+1) 

Piip(n) 

Dp(n  +  1) 

v  (n  +  1) 

pyp(n) 

= 

pyp(n) 

L  r  _ 

y(n  + 1) 

ap(n+  1) 

Equation  (25)  is  significant  because  it  shows  that  “direct  residual  extraction”  is  still  possible  knowing 
only  Qp(n+1)  -  see  equations  (18)  and  (19)  and  associated  discussion.  The  “time  update”  technique, 
presented  above,  foims  an  important  pan  of  the  derivation  of  the  fast  adaptive  filtering  algorithm  pre¬ 
sented  in  this  memorandum.  In  particular,  we  will  explicitly  use  the  decomposition  for  Vp(n)  shown  in 
equation  (26). 
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4.  LINEAR  PREDICTION 


4.1.  Motivation 

The  0(p2)  QRD-based  algorithm  for  the  solution  of  a  p-th  order  adaptive  filtering  problem  devel¬ 
oped  above  has  many  desirable  features  including  being  a  “data  domain”  algorithm  and  having  a  time- 
recursive  formulation,  a  time-independent  computational  requirement  and  a  systolic  architecture  128], 
The  time  shift  redundancy  in  the  adaptive  filtering  problem  can,  however,  be  used  to  reduce  the  com¬ 
putational  load  still  further:  from  0(p2)  to  O(p).  Note  that  the  set  of  rotations,  either  Qp( n)  or  Qp(n), 
that  are  necessary  for  the  solution  of  the  problem  are  entirely  dependent  on  the  matrix  Xp(n).  The  matrix 
X  (n)  can,  however,  be  built  up  in  an  order  recursive  manner  by  adding  extra  columns  which,  because 
oF [he  Toeplitz  nature  of  Xp(n),  consist  of  one  new  element  and  a  time-shifted  version  of  the  previous 
columa  Consider  the  following  decompositions: 


xp(n)  = 


x(l)  ...  x(2 - p) 

x(n)  ...  x(n  -  p  +  1) 


(27) 


where 


=  [xP-i(n>  ^.p-H 

Jxd)  z! 
yf(n)  xp_  j(n  -  1) 


y(<n)  =  [x(2) . x(n))1 


(28) 


(29) 


(30) 


ybp_,(n)  =  (x(2-p) . x(n-p+l)]!  (31) 

and  z  =  [x(0) . x(2-p))'  (32) 

Note  that,  from  equation  (28),  if  we  had  already  determined  the  rotation  matrix  Qp.](n)  that  trian¬ 
gularises  the  matrix6  Xp  j(n)  then  we  could  use  it  to  operate  on  Xp(n)  to  produce  a  partially  triangular- 
ised  matrix.  In  doing  so  we  also  have  to  rotate  the  vector  y^  ^j(n)  however  these  are  exactly  the  steps 

that  are  required  in  the  QRD-based  solution  of  the  (p-l)st  order  backward  linear  prediction  problem.  In 
the  (p-l)st  order  backward  problem,  an  estimate,  at  time  n,  of  x(n-p+l)  is  formed  from  a  linear  combi¬ 
nation  of  the  data  { x(n) . xfn-p+2)).  The  solution  to  this  problem  depends  on  the  triangularisation  of 

the  matrix  Xp  ](n)  and  the  transformation  of  the  reference  vector  y^  pl(n)  -  see  equation  (31).  Hence, 

knowing  the  solution  to  the  (p- 1  )st  order  backward  problem  at  time  n  would  allow  us  to  construct  a  par- 
6.  We  really  mean  the  matrix  B(n)  Xp.|(n)  but  for  the  sake  of  readability  B(n)  will  often  be  omitted. 
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tially  triangularised  version  of  Xp(n)  and  so  save  a  certain  amount  of  computation. 

Equation  (29)  allows  another  partially  triangularised  version  of  Xp(n)  to  be  constructed,  this  time 
using  quantities  from  the  (p-l)st  order  forward  linear  prediction  problem.  The  (p-l)st  order  (forward) 

linear  prediction  problem,  at  time  n,  is  defined  as  the  estimation  of  x(n)  based  upon  the  data  f  x(n-l) . 

x(n-p*l)).  This  involves  the  triangularisation  of  the  matrix  X^j(n-1)  and  the  transformation  of  the  rel¬ 
evant  reference  vector  y^n)  -  see  equation  (30).  As  before  we  can  use  the  decomposition  given  in  equa¬ 
tion  (29)  to  produce  a  partially  triangularised  version  of  Xp(n)  from  known  quantities.  It  should  now  be 
clear  that  the  two  linear  prediction  problems  of  order  (p- 1 )  are  intimately  connected  to  the  problem  of 
determining  a  set  of  rotations  that  triangularise  the  data  matrix  Xp(n).  The  triangularisation  of  Xp(n)  is 
however  central  not  only  to  the  adaptive  filtering  problem  but  also  to  the  p-th  order  linear  prediction 
problems.  We  therefore  have  the  beginnings  Of  an  order  recursive  algorithm  for  linear  prediction  and 
for  adaptive  filtering. 

4.2.  Forward  Linear  Prediction 

The  p-th  order  forward  linear  prediction  problem,  at  time  n,  requires  the  determination  of  the  vec¬ 
tor  of  filter  coefficients  iQfp(n)=  [<  (n),  ....  to'P- '*(n)j  that  minimises  the  total  prediction  error 
Ef  p(n)  =  IIB(n)  £f  p(n)ll  where 


=  Xp(n-l)ffifp(n)  +  y,(n) 


(33) 


As  in  section  3.  the  least  squares  solution  to  this  problem  can  be  found  by  the  method  of  QR  de¬ 
composition.  It  is  necessary  to  determine  the  rotation  matrix  Qp(n-1)  that  triangularises  the  weighted 
data  matrix  B(n-l)Xp(n-l)  and  then  to  apply  it  to  the  weighted  vector  B(n-l)  y^n)  in  order  to  calculate 
the  angle  normalised  residual  af  p(n)  (cf.  equation(19)).  We  also  need  to  be  able  to  calculate  yp(n- 1 ) 
(see  equation(18))  in  order  to  generate  the  a-posteriori  prediction  residual.  Note  also  that  the  triangular¬ 
isation  of  Xp(n- 1 )  is  exactly  what  is  required  in  the  solution  of  the  p-th  order  adaptive  filtering  problem 
at  time  (n-1)  -  see  equation  (4).  Consider,  therefore,  the  following  composite  matrix: 

(V^n)  Xp(n  ~  1)  Jf(n  -  1)  5n  _  ,J  (34) 


From  equations  (22)  and  (23),  it  is  clear  that 

QpOi-l)!^.!*  0p(n-  !)*„_,  =[a‘p(n-l),fi,.Tp(n-l))t  (35) 

where  3p(n-l)  is  a  p-dimensiona)  vector.  It  should  be  clear,  therefore,  that  we  include  the  vector  j  in 
the  above  matrix  (equation  (34))  in  order  to  be  able  to  calculate  yp(n-l)  just  as  the  vector yj(n)  is  present 
to  allow  af  p(n)  to  be  calculated.  Similarly  the  presence  of  the  vector  y(n-l)  will  allow  us  to  calculate 
Op(n-l). 

From  equation  (28)  we  have  that 
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[y/n)  Xp(n  -  1)  y(n  -  I)  sn_ ,]  =  [y/n)  Xp_ ,(n -  1)  yb  p_  ,(n  -  1)  y(n  -  1)  ffji_ ,]  (36) 

Hence  Q^fn-l)  B(n-1)  [yf(n)  Xp_  ,(n  -  1)  yb  p_  ^n-  1)  y(n  -  1)  xn_ ,] 

_  »f. p _  ,(n)  Rp -  i(n- 1)  t(n  - 1)  Bp.  ,(n  -  1 )  3p .  ,(n  -  1  j] 

'[yf.P-i(n)  o  o  v  i("  - 1)  fp. ,("  -  Dj  (37) 

where  Sp.1(n-l)  =  [flt.Yp.1(n-l)]t  (38) 

It  is  clear  that  Vj.p  j(n)  and  p](n-l)  must  have  a  time  recursive  decomposition  similar  to  that 
given  in  equation  (26)  for  Vp  j(n-l).  Hence 


Pf,p-i(n)  Rp_  )(n—  1)  Bb(P_,(n-  1)  Bp_,(n-  1)  Pp.,(n-  lj 

yf.P-i(n)  °  yb.p-i(”-l)yp-i(n-l)f0_.(n-l) 


i(n)  Rp- 

,(n- 

-1)  Bb  p.jfn-  1)  Bp.  j(n  -  1)  9p_  ,(n  —  1) 

(n  —  1 ) 

0 

feb,P-t<n_2)  P^p- t^n  ~2>  <? 

](n) 

o' 

Ob.p-lfo-l)  «p.  j(n-  1)  Yp.,(n-1) 

Now  suppose  that  we  had  already  calculated  a  rotation  matrix,  Qf  p(n- 1 )  say,  that  rotates  the  vector 
yb  p_  j  (n-2)  into  a  form  where  only  the  top  element  is  non-zero.  Then 


O  In  l)ol  Bf'P-'(n)  RP-l("''1)  Pp-,(n-l)  9p.,(n-lj 

f'P  ,  P-vf.p-,(n-D  O  Pyb,p.1<n-2)Pyp.](n-2)  p 

JL  af.P-i(n)  «b.p-,<n-D  op.,(n-l)  Yp_  ,(n  -  1) 

Rp-ifn-l)  Pb,p-t(n-5)  Pp-,(n-l)  flp_ ,(n -  1) 

=  Putp_,<n- 1)  £>'  pEbp_I(n-2)Pnp_1(n-2)  0 

P^-f,  P-  i<n  ~  0  P  P^p_  j(n  -  2)  p 

afp_,(n)  p>  <\p_,(n-l)  ap_,(n- 1)  yp_,(n- 1) 


where  the  new  quantities  Hfj).,(n-1),  ^^(n-l),  iy,(n-2),  2^,  (n-2)  and  ,(n-l)  are  defined  by 
this  operation.  Note  by  analogy  with  equation  ( 1 6)  that  ^  j(iv2)is  the(p-l  )st  order  backward  predic¬ 
tion  energy  at  time  (n-2).  Now  in  order  to  complete  the  iriangularisation  of  the  matrix  X  (n-1)  (see  equa¬ 
tion  (37))  all  that  is  required  is  the  annihilation  of  the  single  element  ot^  ,(n-l).  This  can  be  carried 
out  using  a  single  Givens  rotation: 
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yf,p_,(n)  Rp_ 

i(n  —  1)  J>b>p_,(n-1) 

l>p_i(n-l)  3p_,(n-l) 

Qf,p(n) 

0*  P£b.p-l<n-2> 

Pnp_,(n-2)  0 

PXfp_i(n-l) 

o  0 

Pi-p_,(n-2)  p 

.  atP-i<n> 

“p.^n-D  Yp_i(n-1) 

Uf.p_i(n)  Rp. 

i(n-l)  »b>p_i(n-l)  yp_  j(n  -  1) 

Sp-t(n-l) 

o'  eb,P-i<n-1)  Hp-i(n_1> 

KpOt-D 

P*f.p-i(n-  *) 

O  0  Pip./n-l) 

£> 

1 

JR 

P1  0  ap(n-l) 

Yp(n-l) 

Uf.pfa)  Rp(n  "  1)  yp(n  -  1)  8p(r.  -  1  j 
yfiP(n)  O  yp(n-l)gp(n-l) 


where  Kp(n-l)  is  defined  by  this  operation  and  the  identity  in  equation  (42),  and  hence  the  labelling  of 
some  of  the  elements  in  the  second  matrix  above  follows  by  definition  (see  equation  (37)).  Thus  the  se¬ 
quence  of  orthogonal  transformations  shown  in  equations  (37),  (40)  and  (41)  solve  the  p-th  order  for¬ 
ward  linear  prediction  problem.  Note,  however,  that  the  intermediate  matrix  shown  on  the  right-hand 
side  of  equation  (40)  consists  entirely  of  quantities  that  would  be  available  if  the  (p-l)st  order  forward 
and  backward  problems  had  already  been  solved.  If  this  assumption  were  true  then  we  could  have  con¬ 
structed  this  intermediate  matrix  directly,  thereby  circumventing  the  need  for  the  operations  as  outlined 
in  equations  (37)  and  (40).  Only  the  single  Givens  rotation  of  equation  (41 )  need  actually  be  performed 
This  requires  a  fixed  amount  of  computation  because  only  eight  elements,  one  of  which  is  zero,  of  the 
left  hand  matrix  in  equation  (41)  are  affected  by  the  required  rotations.  In  order  to  complete  the  lattice 
algorithm  for  the  linear  prediction  problem,  the  fast  update  for  the  auxiliary  (backward)  problem  must 
be  derived.  This  can  be  done  along  similar  line  to  the  above  (sec  section  4.3). 

Before  considering  the  backward  linear  prediction  problem,  note  from  above  that  the  "new”  quan 
tities  Kp(n-l)  and  i^p.j(n-l)  have  the  following  interpretation: 

^f,p-l(n'1)  =  *f,p(n'1)  <4?) 


V,(n-1) 
_  Vn_  1). 


=  2p(n-l) 


(44) 


Note  also  that  equation  (41)  provides  a  recursive  decomposition  of  the  matrix  Rp(n-I).  This  show  s 
that  the  diagonal  elements  of  this  matrix  are  in  fact  the  backward  prediction  residual  energy  terms  for 
each  of  the  sub-order  problems.  It  also  explains  why  the  Givens  rotations  used  in  QRD-based  linear  pre¬ 
diction  algorithms  should  ensure  that  the  diagonal  elements  of  the  R  matrix  arc  positive  (see  Bellangcr 
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and  Regalia[20]  for  further  insight). 

4.3.  Backward  Linear  Prediction 


The  p-th  order  backward  linear  prediction  problem,  at  time  n,  requires  the  determination  of  the  vec¬ 
tor  of  filter  coefficients  ^(n)  =  [co^fn),  ....  CD^pp~  !)(n)j  that  minimises  the  total  prediction  error 
Eb  pOO  =  »B(n)  £^(n)ll  where 

fib.p(")  =  Xp<n)  +  ib,p<n>  <45> 


Again  the  least  squares  solution  to  this  problem  can  be  found  by  the  method  of  QR  decomposition. 
It  is  necessary  to  determine  the  rotation  matrix  Qp(n)  that  triangularises  the  data  matrix  Xp(n)  and  then 
to  apply  it  to  the  vector  yb  p(n)  in  order  to  calculate  p(n)  (cf.  equation  (19».  We  also  need  to  be  able 
to  calculate  yp(n)  (see  equation  (1 8))  in  order  to  generate  the  a-posteriori  prediction  residual.  Consider, 
therefore,  the  following  composite  matrix  and  the  illustrated  decomposition: 


O1  0  0 

l(n_1)l'b.p-i(n-1)5n-l_ 


(46) 


In  equation  (46),  it  has  been  assumed  that  the  data  sequence  x(n)  is  pre-windowed  (i.e.  x(n)  =  0  for 
n  5  0).  Note  that  this  is  the  only  place  in  the  analysis  that  requires  this  assumption.  Consider  the  effect 
of  the  rotation  matrix  Qp _1(n-l)  on  the  lower  n-1  rows  of  the  matnx  in  equation  (46)  -  after  weighting 
by  B(n)  of  course: 


1  °l  B(n)  0  0  0 

P  Qp-  |("  -  »J  I?/")  Xp-  l(n~  !)  ?b.p-  1(n  "  l)  *n-  1 

P"  ~ 1  x(  1 )  o'  0  0 

=  i'r,p-i(n)Rp_1(n-l)ybiP.)(n-l)sp.1(n-l)  (47) 

p-  l(n)  0  yb.p-t^-Dfp. !("-») 


As  before,  all  the  vectors  on  the  bottom  row  of  the  above  matrix  have  a  decomposition  in  terms  of 
their  value  at  the  previous  time  instant  and  a  new  element: 
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P*"‘x(l) 

p1 

0 

0 

!(n)  Rp_ 

l<n“ 

Dw^p.idt- 

l)sp- 

(n-1) 

yf.P-i(n) 

0 

yb.P-i("- 

l)h- 

(n-1) 

“-*x(l) 

o' 

0 

0 

“ 

RP 

.,(n- 

•n  »b,P-i(n 

-1)  *p 

-|(n" 

1) 

Py^-D 

O 

P^b.p-i<n 

-2) 

P 

L  « 

CP-l(n> 

0‘ 

“b.p-t(n 

-1)  Yp 

1) 

(48) 


Now  suppose  that  we  have  already  constructed  a  rotation  matrix  QbD(n-l)  that  annihilated  the 

n-2  *P 


vector  yf  ^j(n-l)  by  rotation  against  the  element  pn’^x(l).  Then 


Qb  p(n-  l)o 


^-’xd) 

o'  0 

«t,p-t<n)  RP- 

,(n-1)  pb  p_j(n-  1)  a, 

°  Pyb,p-i(n-2) 

af.P-tW 

«b,p-,(n-1)  Y, 

p' 

ppb.p.1(n-2)  0 

l»tp-,(n)  Rp. 

,(n- 

D  "b.p-iC1"1)  9p-,(n-l) 

0 

O 

P^b,P-](n'2)  Q 

.  af,P-i<n> 

P1 

“b.p-t^"1)  Yp-,(n-  1) 

(49) 


Let  Qb  p(n)  be  the  rotation  matrix  that  annihilates  the  element  a(  ^(n)  by  rotation  against  the  el¬ 
ement  fkf  p  ](n- 1 ).  Application  of  the  transformation  Qb  p(n)  to  the  above  matrix  thus  yields  the  result: 


Qb.p(nJ 


Petr-i<n  ~ 9'  Ph,,p-i<n~2>  0 

P  O  P^b,p-i(n_2)  P 

i(n)  o' 


“b.p-)^-1)  Yp.,(n-1) 


hv-W  0  ^b.P-i<n- *) 


<Pp(n) 


Utp_,(n)Rp.,{n-l)  l>b  p.](n-  1)  Sp.^n-l)1 


0 

0 


O  P^b.p_,(n-2) 


b,  p' 


(50) 
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where  the  new  quantity  <Pp(n)  is  defined  by  this  equation.  The  identity  of  the  elements  in  the  bottom  row 
of  the  above  matrix  follows  because  of  the  following  reasoning:  bearing  in  mind  the  underlying  data 
matrix  (see  equation  (46)),  we  are  attempting  to  create  an  upper-triangular  pxp  matrix  in  the  upper  left- 
hand  comer  on  the  matrix  in  equation  (SO).  At  present  this  sub-matrix  is  rather  sparse  and  a  little  thought 
shows  that  it  is  easy  to  construct  a  matrix  (Qb  p(n)  say  )that  will  perform  this  triangularisation7.  Spe¬ 
cifically. 


o' 

Hb.p- 

!<"- 

1) 

9p(n) 

Qb,  p(n) 

Rp-ifa- 

*>  ub.p- 

l("- 

»  8p 

_,(n-  1) 

P 

o 

p*b.r- 

|(n- 

-2) 

P 

0 

Dl 

«b 

p(n) 

Vn) 

Rp(n)  ubi 

p(n)  V") 

•O 

>1 

o 

p(n)  gp(n) 

The  important  thing  to  note  is  that  the  action  of  the  matrix  Qb  p(n)  on  the  above  matrix  only  af¬ 
fects  the  upper  p  rows  of  elements.  Thus  we  conclude  that  the  lower  n-p  rows  of  elements  of  the  matrix 
in  equation  (50)  must  be  equal  to  the  lower  n-p  rows  of  elements  of  the  matrix  in  equation  (5 1 )  and  the 
result  holds. 

Hence  the  sequence  of  orthogonal  rotations  given  in  equations  (47),  (49),  (50)  and  (51)  solve  the 
p-th  order  backward  linear  prediction  problem.  Following  the  development  of  the  solution  to  the  for¬ 
ward  problem  in  section  4.2,  note  that  the  data  matrix  on  the  left-hand  side  of  equation  (50)  can  be  con¬ 
structed  directly  given  the  solutions  to  the  (p- 1  )st  order  forward  and  backward  problems.  Thus  the  trans¬ 
formations  shown  in  equations  (47)  and  (49)  can  be  by-passed.  Furthermore,  as  both  p(n)  and  yp(n) 
are  available  in  the  matrix  on  the  right-hand  side  of  equation  (50),  the  transformation  shown  in  equation 
(51)  is  not  required  either,  provided  we  are  only  interested  in  the  prediction  residuals.  Thus  only  the  ro¬ 
tations  summarised  in  equation  (50)  need  actually  be  performed  explicitly.  This  requires  a  fixed  amount 
of  compulation  because  only  six  elements,  one  of  which  is  zero,  of  the  left  hand  matrix  in  equation  (50) 
are  affected  by  the  required  rotations. 


7.  In  actual  fact,  it  is  not  necessary  that  we  specify  how  this  triangularisation  is  achieved  since  the  QR 
decomposition  theorem  guarantees  us  that  it  is  possible.  Note,  however,  that  it  is  crucial  to  QRD-bascd 
fast  Kalman  algorithm  that  this  transformation  is  done  in  a  specific  way.  Interestingly,  it  has  recently 

been  pointed  out  [20]  that  the  sines  of  the  angles  involved  in  the  Qbi  p(n)  rotation  are  in  fact  the  re¬ 
flection  coefficients  of  a  "conventional''  least  squares  lattice  algorithm. 


16 


Note,  in  passing,  that  the  "new"  quantities  tp^fn)  and  ^(n-2)  have  the  tollowtng  interpretation. 

Vi(n'2)=Vn,)  <52> 


<P  (n) 


Gathering  lOgether  the  results  of  sections  4.2  and  4.3  we  see  that  it  is  possible  to  transform  various 
quantities  from  the  solution  to  the  (p-l)st  order  forward  and  backward  linear  prediction  problems,  at 
time  n  and  (n-1)  respectively,  into  the  same  quantities  from  the  solution  to  the  p-th  order  problems  at 
time  n  (see  Figure  2).  Thus,  given  that  0-th  order  linear  prediction  is  trivial,  we  can  generate  the  solution 


to  the  p-th  order  problem  by  iteration  in  order.  The  resultant  architecture  has  a  lattice  structure  and,  since 
the  number  of  operations  per  stage  is  fixed,  requires  0(p)  operations.  Note  that  by  including  the  adaptive 
filtering  reference  vector  y  (n-1)  in  the  calculation  of  the  p-th  order  forward  linear  prediction  problem 
(section  4.2)  we  automatically  solve  the  p-th  order  adaptive  filtering  problem  for  time  (n-1 ). 

In  most  cases,  the  fact  that  the  adaptive  filtering  residual  is  delayed  by  one  time-step  will  be  of  no 
consequence,  especially  given  the  regularity  of  the  algorithm/architecture  and  the  benefits  this  brings 
with  it.  It  is  possible,  however,  to  alter  the  algorithm  presented  above  so  as  to  solve  the  adaptive  filtering 
problem  at  time  n  but  at  the  expense  of  more  complicated  mathematics  and  a  slightly  less  regular  archi¬ 
tecture:  see  section  10.1  for  more  details. 

Note  from  equation  (41)  and  equations  (50)  and  (51)  that  it  is  possible  to  transform  a  matrix  of 
quantities  from  the  solution  to  the  (p- 1  )st  order  forward  problem  at  time  n  and  a  similar  matrix  from  the 
solution  to  the  (p-l)st  order  backward  problem  at  time  n  into  the  same  quantity  (viz  Rp(n)).  Thus  it  is 
possible  to  transform  quantities  from  the  forward  problem  at  time  n  into  quantities  from  the  backward 
problem  at  time  n,  although  this  involves  a  numerically  unsound  “inverse  rotation".  This  is  one  of  the 
crucial  steps  in  the  derivation  of  the  QRD  fast  Kalman  algorithml3]  (see  Figure  3).  The  remaining  step 
is  to  deduce  the  rotation  Qp(n+1)  from  knowledge  of  quantities  from  the  backward  problems  at  time  n 
and  (n-1).  This  latter  step  involves  using  the  pinning  vector  as  in  the  classical  Fast  Kalman  algorithm 
(see  Proudler,  Me  Whiner  and  Shepherd!  16]  for  more  details). 
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5.  IMPLEMENTATION 


The  QRD-based  lattice  algorithm  developed  above  can  be  implemented  using  the  processor  archi¬ 
tecture  shown  in  Figure  4.  The  functionality  of  the  processing  elements  (PE’s)  is  shown  in  Figure  S  and 
a  “pseudo-code”  listing  of  this  algorithm  is  given  in  section  10.3. 

In  the  derivation  of  the  algorithm  (section  4),  no  mention  was  made  of  how  the  Givens  rotations 
were  to  be  implemented.  There  are  several  ways  in  which  this  can  be  done.  The  obvious,  or  “normal”, 
approach  involves  the  calculation  of  a  square-root  (see  Figure  3),  which  is  a  computationally  intensive 
step  when  using  floating  point  numbers.  It  is  possible  to  implement  the  square  root  operation  efficiently 
by  use  of  CORDIC  algorithms[27].  However  this  requires  the  use  of  fixed  point  arithmetic  and  it  has 
been  found  [28]  that  for  most  practical  applications,  a  floating  point  implementation  is  required.  On  the 
other  hand,  Gentleman[6]  and  Hammariing[7]  have  derived  a  modified  Givens  rotation  that  requires  no 
square-root  operation  and  this  is  clearly  advantageous  for  floating  point  implementations. 

The  essence  of  this  “square-root-free”  technique  is  to  factorise  the  triangularised  matrix  (e.g. 
Rp(n»  into  two  parts,  one  of  which  contains  ail  the  quantities  which  involve  a  square-root.  This  laner 
factor  can  be  calculated  indirectly  as  its  square  thus  avoiding  the  square-root  operation  (see  section 
10.2).  However  this  means  that  the  rotation  (an  orthogonal  transformation)  is  now  implemented  by  two 
separate  non-orthogonal  transformations  acting  upon  the  factorised  quantities.  As  such,  one  would  be 
quite  justified  in  questioning  whether  or  not  this  square-root-free  Givens  rotation  leads  to  an  orthogonal 
algorithm.  Nevertheless,  simulations  have  shown  that,  in  floating  point  arithmetic  at  least,  the  squarc- 
root-free  version  is  more  stable  numerically  (see  section  7).  This  is  despite  the  reservations  about  the 
“square-root-free”  Givens  rotations  algorithm  that  have  been  raised  by  the  numerical  analysis  commu¬ 
nity  [6][7], 

Closer  analysis  of  the  normal  Givens  rotation  shows  that  square-root  operation  is  required  only  in 
situations  where  the  square-root  of  the  sum  of  the  squares  of  two  quantities  is  required.  The  numerical 
stability  of  the  “square-root-free”  method  may  well  be  due8  to  the  fact  that  the  act  of  squaring  these  two 
numbers  only  to  take  the  square-root  of  their  sum  is  numerically  inferior  to  propagating  the  squared 
quantities  directly. 

The  “square-root-free”  algorithm  mentioned  above  is  an  algorithm  for  calculating  the  required  pla¬ 
nar  rotation.  In  the  QRD-based  least  squares  minimisation  problem  these  rotations  also  have  to  be  ap¬ 
plied  to  various  vectors.  It  is  possible  to  implement  a  rotation  in  two  ways:  in  a  feedforward  or  a  feed¬ 
back  mode.  When  implementing  a  planar  rotation  (a  two-input,  two-output  transformation),  the  normal 
choice  would  be  to  calculate  each  output  in  terms  of  the  two  inputs  separately.  This  leads  to  a  “feedfor¬ 
ward"  system.  However,  it  is  possible  to  reformulate  the  transformation  so  that  one  output  is  now  de¬ 
pendent  on  one  input  and  the  other  output  Initially  the  motivation  behind  this  reformulation  was  to  re¬ 
duce  the  number  of  multiplications  required  [23].  It  has  been  shown,  however,  that  this  alternative  im¬ 
plementation  corresponded  to  the  "error-feedback”  technique  proposed  by  Ling  and  Proakis  [12]. 

It  is  believed  that  this  feedback  has  a  stabilising  effect  when  errors  are  made  in  the  calculation  due 
to  finite-precision  effects  in  the  arithmetic.  Indeed,  computer  simulations  (see  section  7)  have  shown 

8.  P  A  Regalia,  Private  Communication. 
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that  this  error-feedback  has  a  significant  effect  on  the  numerical  stability  of  the  QRD-based  lattice  al¬ 
gorithm:  the  feedback  version  is  quite  capable  of  working  with  4  bit  mantissas  (at  least  for  sequence 
lengths  of  up  to  10,000  -  the  longest  simulation  tun  so  far).  It  is  interesting  to  note  that  the  feedback 
structure  can  be  thought  of  as  one  half  planar  rotation  and  one  half  hyperbolic  rotation.  It  is  well  known 
that  a  planar  rotation  is  numerically  superior  to  a  hyperbolic  rotation;  however,  it  would  appear  that  this 
half-and-half  mixture  is  better  than  both. 

Combined  with  the  two  methods  for  calculating  the  rotation  parameters,  the  feedforward/feedback 
choice  results  in  four  different  variations.  The  computer  simulations  of  section  7  show  that,  of  these  four 
possible  variations,  the  one  that  is  equivalent  to  the  RMGS  error-feedback  algorithm  (square-root-free 
with  feedback)  performs  the  best  in  terms  of  numerical  stability.  A  “pseudo-code”  listing  of  this  version 
of  the  algorithm  is  given  in  the  appendix.  It  is  worth  emphasising  that  the  basic  architecture  (Figure  4) 
is  not  affected  by  the  choice  of  rotation  technique  so  that  the  only  difference  between  the  different  op¬ 
tions  is  a  change  of  PE's.  The  function  of  the  "square-root-free  with  feedback”  PE’s  is  shown  in  Figure 
6.  Comparison  with  the  Ofp2)  QRD-based  algorithm  [28]  shows  that  the  PE’s  presented  here  are  exactly 
the  same  as  used  in  the  triangular  systolic  array.  Indeed,  when  the  QRD-based  lattice  algorithm  is  gen¬ 
eralised  to  the  multi-channel  case  (section  6),  the  processing  required  in  each  of  the  lattice  stages  neces¬ 
sitates  the  use  of  triangular  arrays  of  PE’s.  In  fact  the  lattice  stages  shown  in  Figure  4  may  be  considered 
to  consist  of  lxl  triangular  arrays! 
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Delayed  Adaptive 
Filtering  Residual 


Forward 

Residual 


Backward 

Residual 


"T*  -*0-— 

Multiplier  Delay 

Figure  4‘  Delayed"  Adaptive  Filtering  Lattice 


e'  :=J(KP)2  +  \ayJ2- 

IFe'  =  0  THEN  cx  p+1  ~Lsx  p+1  :=0 
ELSE  cx  p+j  •=  P^  p  !  £  •  s*,p+]  =  “y.p  !  c 
END_IF; 

Yp+1  •“  cx,p+l  Yp- 

l1  ,=  cx,p+l  P^X.p  +  sx.  p*  1  ^.p' 

“x.p+l  ;=  cx.p+l  “x.p  •  SX,p+l  P^x.p; 


Figure  5  "Normal"  Rotation  PE's 
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6.  MULTI-CHANNEL  CASE 


6.1.  Problem  Definition 

The  extension  of  the  adaptive  filtering  algorithm  presented  above  to  the  wide-band  bcamforming 
problem  is  relatively  straight  forward:  all  that  is  require  is  that  certain  scalar  quantities  be  replaced  by 
vectors  and  some  vectors  be  replace  by  matrices.  The  essential  features  of  the  derivation  presented  in 
section  4  carry  over  exactly.  Rather  than  reproduce  the  mathematics  of  section  4,  in  the  following  we 
merely  outline  the  salient  points  of  the  derivation  of  the  multi-channel  algorithm.  As  before,  we  consid¬ 
er  only  real  signals:  the  extension  to  the  complex  case  is  straight  forward. 

In  a  multi-channel  least  squares  adaptive  filtering  problem,  a  set  of  N  p-dimensional  weight  sec¬ 
tors,  Jpp'*(n)  £0  <  i  <  N-l),  is  to  be  found,  at  time  n,  that  minimises  that  the  sum  of  the  differences, 
squared,  between  a  reference  signal  y(n)  and  a  linear  combination  of  N  samples  from  each  of  p  data  time 
series  x/n-j)  (1  S  i  S  p,  0  S  j  S  N-l).  This  is  equivalent  to  adaptively  filtering  p  separate  time  series  in 
order  to  form  the  best  estimate  of  the  reference  signal  (see  Figure  7).  If  the  p  data  sequences  come  from 
spatially  separate  antennae  then  we  have  a  spatial  as  well  as  a  temporal  filtering  problem.  Specifically . 
the  measure  IIB(n)  gN(n)ll  is  to  be  minimised,  where: 

£N(n)  =  XN(n)aN(n)  +  y(n)  (54) 
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XN(n)  = 


X,(l)  X2{1)  .. 

•  V1) 

«[(0)  Xj(0)  . 

-  «p(0) 

x,(2-N)  x2(2-N)  ...  xp(2  -  N) 

*j(n)  x2(n)  ., 

.  *p(n) 

x,(n-l)  x2(n- 1)  .. 

■  Xp(n-l)  . 

Xj(n  -  N  +  1)  x2(n  -  N  +  1)  ...  xp(n  -  N  ♦  1) 

(55) 


x‘(l)...  x'(2  -  N) 
»‘(n)  ...  X‘(n-N+1) 


(56) 


fi>N(n)  = 


iPp0)(n) 


JPpN-1)(ny 


(57) 


and 


y(n)*fy(l). ....  yfn)]1 


(58) 


Equation  (56)  serves  to  define  the  new  vector  quantity  x(n).  Note  that,  apart  from  the  change  from 
scalar  to  vector  quantities,  equation  (56)  is  identical  to  equation  (5).  The  solution  of  this  least  squares 
minimisation  problem  via  QR  decomposition  is  no  different  from  any  other:  it  requires  the  determina¬ 
tion  of  an  orthogonal  matrix  QN(n)  that  transforms  the  matrix  B(n)  XN(n)  into  an  uppenriangular  form 
The  fact  that  the  matrix  XN(n)  is  block-Toeplitz  allows  us  to  use  the  ideas  developed  in  section  4  to 
construct  an  order  recursive  algorithm.  Again  this  relies  on  the  iterative  structure  present  in  the  multi¬ 
channel  linear  prediction  problems. 

6.2.  Forward  Linear  Prediction 

In  the  N-th  order  multi-channel  forward  linear  prediction  problem,  an  estimate  of  the  vector x(n) 

is  formed,  at  time  n,  from  a  linear  combination  of  the  data  (x(n-l) . x(n-N))  .  Thus  it  is  necessary  to 

determine  the  rotation  matrix  QN(n-l)  that  triangularises  the  data  matrix  XN(n-l).  Consider,  therefore, 
the  following  composite  matrix: 


[Y,(n)  XN(n-  1)  y(n  -  1) 


(59) 


where  the  "reference  matrix"  Yf(n)  is  defined  as 
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Yj(n)  s 


x'(2) 

x,(2)  x2(2)  . 

•  xp(2) 

An) 

Xj(n)  x2(n)  . 

•  xp(n) 

(60) 


and  it  the  multi-channel  equivalent  to  if(n)  -  see  equation  section  (30).  From  equation  (56)  we  have  that 
[Yj(n)  XN(n-  1)  y(n-  1)  =  [y,( n)  XN_,(n-  1)  Yb >N_,(n-  1)  y(n-  1)  (61) 


where  Y^.iCn-l)  is  the  reference  matrix  for  the  (N-l)st  order  backward  linear  prediction  problem  at 
time  (n-1)  (see  section  4.3): 


x‘(2  -  N) 

= 

x,(2  -  N)  x2(2-N)  . 

-  xp(2-N) 

An  -  N) 

Xj(n  -  N)  x2(n  -  N)  . 

.  xp(n-N) 

Hence  QN_,(n-l)  B(n-l)  [Y^n)  XN.,(n-  1)  Yb  N_  ,(n-  I)  y(n-  1)  Sn_,] 

(ff, N _  j(n)  Rjm _  i(n  —  1)  Ub  n_  j(n  -  1)  pN_  j(n  -  1 )  a*,- _  ](n  —  1)| 
'  vfiN-i(n)  O  Yb  n_  j(n-  1)  yN.  ,(n  -  1)  fN  ,(n  -  l)j 


where  the  matrices  Ufjs'  i(n),  Vf  N.,(n),  UbN.|(n-l)  and  VbN  j(n-l)  arc  the  multi-channel  equivalents 
of  Uf.p-i(n),  Vf  N.i(n),  Uh.p  i(p-l)  and  Ybp  l(p- 1)  respectively.  Note  that  Vf  N.,(n)  and  VbN.,(n-l)  have 
a  time  recursive  decomposition  similar  to  that  given  in  equation  (26)  forv^  jtn-1)  and  that  RN.>(n-l)  is 
a  (N-l)px(N-l)p  upper  triangular  matrix. 

Now  suppose  that  we  had  already  calculated  a  rotation  matrix,  Q^Oi-l)  say,  that  transforms  the 
matrix  j  (n-2)  into  an  upper-triangular  form9.  Then 


Qf,N(n-  1)  o 


Rn-i(d  —  ^b,N-I^0  — ^  N  —  1  ^ 

PVf,N_  i(n  —  1)  O  pVfc  N_  j(n  -2)  PyN  _  j(n  -  2)  P 

ff'f,N_i(n)  f>‘  »h,N-](n“  D  >)  »_ 


pMfiN_,(n-l) 
PAf ,N_j(n-l) 

fff,N-l(n) 


O 

o 

-l 


PEb,N-l<D_2)  PyN-l(n-2)  P 
o  pXN.,(n-2)  p 

«U-l(r,_1)  YN_,(n-l) 


(64) 


9.  Remember  that  we  are  restricted  to  using  only  orthogonal  operations! 
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Note  that  there  are  now  p  residuals  for  both  the  forward  and  backward  linear  prediction  problems  hence 
the  vectors  fflf,N.i(n)  and  SbjM.  t(n-l ).  Now  in  order  to  complete  the  triangularisation  of  the  matrix  XN<n- 
1)  all  that  is  required  is  the  annihilation  of  the  vector  g^iCn-l).  This  can  be  carried  out  using  a  se¬ 
quence  of  Givens  rotations  rather  like  that  used  in  equation  (23): 


pMfN_,(n-l) 

PAfN_,(n-l) 


Uf,N-lOO  RN- 


M 


f,N 


,i(n) 
PAf,N_,(n-D 

»f.N(n) 


Rn-i^**— J)  Un _ i(n  8)>j_](n—  1) 

°  PEb.N-lfB-2)PilS-l(n-2)  0 

O  O  pX.N_,(n-2)  o 

P1  fi^N-t^-1)  aN-t(n_1>  rN-i<n-i) 

j(h  O^N-t^"!)  1)  d|s4_|(ri  1) 

0  Eb,  n  -  l(n  ~  KN(n-l) 

O  O  pXN_  ,(n  -  2)  9 
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Note  that  the  intermediate  matrix  shown  on  the  left-hand  side  of  equation  (41)  consists  entirely  of 
quantities  that  would  be  available  if  the  (N-l  )st  order  forward  and  backward  problems  had  already  been 
solved.  Thus  only  the  Givens  rotations  defined  by  equation  (41)  need  actually  be  performed.  This  can 
be  implemented  on  the  standard  triangular  systolic  array  [28]  in  0(p2)  operations  per  time  instant. 

6.3.  Backward  Linear  Prediction 

The  N-th  order  backward  linear  prediction  problem  is  defined  as  the  estimation,  at  time  n.  of 

x(n-N)  from  a  linear  combination  of  the  data  (x(n) . x(n-N+l)}.  Consider,  therefore,  the  following 

composite  matrix  and  the  illustrated  decomposition: 


[Xp(n)  Yb,N(n)  5b]  = 
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The  rotation  matrix  QVI(n-l)  will  triangularise  the  data  matrix  XN  I(n-l)  hence: 
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where  we  have  explicitly  used  the  familiar  decomposition  for  the  matrices  Vfjj.,(n)  and  Vb  N.,(n-l). 

The  calculation  of  the  quantities  necessary  for  direct  residual  extraction  can  now  be  achieved  in 
three  steps.  First  imagine  that  we  had  operated  on  the  above  matrix  with  an  orthogonal  transformation 
(Q<>  -  1))  that  moves  the  top  row  of  elements  down  to  the  row  just  above  the  top  of  the  matrix 
vr.N-i(n-l)-  The"  suppose  that  we  have  already  constructed  a  rotation  matrix  Q^-fn-l)  that  performed 
a  QR  decomposition  on  the  composite  matrix 
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Finally,  following  section  4.3,  let  Qb_\(n)  be  the  rotation  matrix  that  annihilates  the  row  vector 
C’f,  N  .  j(n)  by  rotation  against  the  triangular  matrix  (3Ef  N1(n- 1 )  Then 


Uf,yt_  j(n) 

Rn 

_t(n 

-1) 

Ub,N-t 

(n- 

-»  »N- 

,(n- 

1) 

/ft  Y 

PEf.N- |(n“ 

1) 

0 

PMb.N- 

t(n 

-2) 

P 

O 

0 

PAb,  N  - 

,(n 

-2) 

P 

«IN- 1<"> 

P1 

®U-l 

(n- 

-1)  Yn- 

,(n- 

1) 

u 

f,N_i(n)  Rn_ 

,(n- 

1)  u, 

b,  N  - 

,(n  -  1) 

*N 

_  ,(n  -  1)' 

E 

r,N-i<n) 

0 

M 

b,  N  - 

.  j(n  —  1 ) 

S*N(n) 

0 

0 

P  A 

‘b,N 

_  ,(n  -  2) 

P 

P1 

P’ 

ffb, 

n("> 

YN(n) 

26 


Delay 


Triangular 

Swith 
a 


Figure  8.  Lattice  of  Triangular  Arrays. 


and  we  have  achieved  our  objective.  As  in  section  4.3,  we  do  not  have  to  complete  the  lormation  ol  the 
triangular  matrix  (an  NpxNp  matrix  in  this  case)  since  any  operation  necessary  to  achieve  this  would 
not  affect  the  lower  two  rows  of  the  composite  matrix  shown  in  equation  (50).  Again  the  rotation  shown 
in  equation  (50)  can  be  implemented  using  a  triangular  systolic  array  in  0(p2)  operations 

Thus  by  the  use  of  equations  (41)  and  (50),  we  can  calculate  the  solution  to  the  N-th  order  multi¬ 
channel  linear  prediction  problems  in  Otp2)  operations  given  the  solutions  to  the  (N-l)st  order  prob¬ 
lems.  As  before,  the  rotation  that  are  necessary  to  solve  the  forward  linear  prediction  problem  automat¬ 
ically  solve  the  p-ih  order  adaptive  Tillering  problem  for  time  (n-1)  as  well.  The  resultant  architecture 
(see  Figure  8)  has  a  lattice  structure  where  each  stage  of  the  lattice  contains  two  triangular  systolic  ar¬ 
rays  (Figure  9).  Each  of  these  arrays  solve  a  p-th  order  recursive  least  squares  minimisation.  The  total 
number  of  operations  necessary  to  solve  an  N-th  order  multi-channel  adaptive  filtering  problem,  with  p 
channels,  is  0(Np2). 


Once  again,  it  is  possible  to  alter  the  algorithm  presented  above  so  as  to  solve  the  adaptive  filtering 
problem  at  time  n  but  at  the  expense  of  more  complicated  mathematics  and  a  slightly  less  regular  archi- 


tenure:  see  section  10.1  for  details  of  the  single  channel  case:  the  extension  to  the  multi-channel  case 
should  be  obvious. 


7.  SIMULATIONS 


The  single-channel  QRD-based  lattice  filter  algorithm  has  been  simulated  on  a  computer  and  com¬ 
pared  with  the  full  0(p2)  QRD  triangular  array  algorithm.  The  results  show  that  the  two  algorithms  pro¬ 
duce  virtually  identical  results  for  "sensible"  wordlengths.  As  the  wordlength  is  reduced  the  lattice  al¬ 
gorithm  begins  to  suffer  before  the  full  triangular  array  does.  This  is  to  be  expected  as  the  lattice  algo 
rithm  relies  on  an  exact  mathematical  relationship  between  the  forward  and  backward  linear  prediction 
problems  which  is  progressively  made  void  as  numerical  errors  are  made. 

The  situation  considered10  is  that  of  a  channel  equaliser  (Figure  10)  for  a  data  channel  that  has  a 
“raised  cosine”  impulse  response  (equation  (71)). 


n  =  1, 2,  3. 
otherwise 


(71) 


By  varying  the  parameter  W,  the  amount  of  inters ymboi-interfenence  between  a  given  symbol  and  the 
two  either  side  of  it  can  be  changed.  This  in  effect  controls  the  eigenvalue  spread  of  the  data  covariance 
matrix  (see  Table  1 .). 

An  1 1th  order  adaptive  QRD-based  least-squares  lattice  filter  is  used  to  equalise  the  channel  re¬ 
sponse.  In  order  to  “train"  the  equaliser,  the  transmission  channel  is  fed  with  a  polar  (±1 )  pseudoran¬ 
dom  training  sequence.  This  sequence,  delayed  by  seven  time  instants,  is  used  as  the  reference  signal 
for  the  adaptive  filtering  algorithm .  The  delay  is  inserted  to  ensure  that  the  adaptive  filter  has  an  impulse 

10.  Suggested  by  S  Haykin. 
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response  that  is  symmetrical  about  the  centre  tap.  A  small  quantity  of  “measurement'’  noise,  in  the  form 
of  a  pseudorandom  sequence  with  an  approximately  gaussian  probability  distribution  function,  is  added 
to  the  channel  output  The  noise  sequence  used  has  zero  mean  and  a  variance  of  0.001 .  The  “forget  fac¬ 
tor"  p  (see  equation  (2))  was  fixed  at  a  value  of  0.996,  which  implies  an  effective  data  window  of  250 
time  samples. 

All  calculations  within  the  algorithm  were  performed  using  limited-precision  floating  point  arith¬ 
metic.  Only  the  number  of  bits  in  the  mantissa  are  varied  in  the  experiments:  the  number  of  bits  in  the 
exponent  is  fixed  at  eight.  No  quantities  internal  to  the  adaptive  filtering  algorithm  are  held  to  a  greater 
precision  than  for  the  outputs:  the  results  of  all  arithmetic  operations  are  immediately  reduced  to  the 
required  precision. 

The  performance  of  the  equaliser  is  monitored  by  recording  the  ensemble-averaged,  squared  a-pri¬ 
ori  equalisation  error  (see  equation  (10)).  This  has  the  advantage  that  it  shows  how  close  to  convergence 
the  algorithm  is  whilst  still  showing,  asymptotically,  the  least  square  equalisation  error  The  ensemble 
average  is  taken  over  100  realisations  of  the  experiment.  Several  experiments  were  performed  using 
various  combinations  of  parameters  and  algorithms  The  main  results  are  discussed  below,  however,  for 
completeness  a  complete  set  of  results  are  reproduced  in  the  attached  annex. 

Figure  1 1  and  Figure  12  show  the  basic  performance  of  the  QRD-based  lattice  equaliser  system  for 
different  values  of  wordlength  and  eigenvalue  spread.  Figure  11  shows  that,  with  double-precision 
arithmetic,  the  rate  of  convergence  is  more  or  less  insensitive  to  the  different  eigenvalue  spread  settings: 
as  would  be  expected  from  a  recursive  least  squares  minimisation  process.  Figure  12  illustrates  how  the 
wordlength  affects  the  performance  for  a  fixed  eigenvalue  spread  (W=2.9).  Note  that  there  is  very  little 
discemable  difference  between  the  systems  using  12, 16  and  56  (IEEE  double-precision)  bit  mantissas. 

Figure  13  shows  a  comparison  of  the  QRD-based  lattice  algorithm  with  the  full  QRD-based  trian¬ 
gular  systolic  array  version.  Two  other  systems  are  also  shown  in  this  figure:  they  are  the  “square-root- 
free  with  feedback"  forms  of  the  lattice  algorithm  and  the  array  algorithm.  This  figure  shows  the  case 
of  A  bit  mantissas  and  a  fixed  eigenvalue  spread  setting  (W=2.9).  This  may  be  considered  to  be  an  ex¬ 
cessively  short  wordlength.  The  reason  for  this  choice  is  that  often  finite  precision  effect  are  often  only 
manifested  after  the  round-off  errors  have  had  time  to  accumulate  [2],  By  using  a  small  wordlength.  the 
appearance  of  such  effects  occur  sooner  thus  reducing  the  time  necessary  to  perform  the  simulation. 

In  most  cases,  a  p-th  order  adaptive  filter  will  converge  within  2p  time  instants.  At  this  point  the 
a-priori  residual  will  have  reached  a  value  primarily  determined  by  the  eigenvalue  spread  and  not  the 
wordlength.  As  round-off  errors  accumulate,  the  a-priori  error  will  increase  indicating  a  loss  of  accuracy 
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in  the  algorithm.  Up  to  a  run  length  of  UXXX),  the  longest  simulation  run  to  date,  the  normal  lattice  al¬ 
gorithm  retains  its  post-convergence  accuracy  for  12  bit  mantissas.  In  the  case  of  the  square-root-free 
with  feedback  lattice,  the  same  behaviour  is  seen  using  only  4  bit  mantissas. 


It  can  be  seen,  in  Figure  1 3,  that  the  faster,  lattice  algorithm  is  only  marginally  worse  than  the  full 
triangular  array  version  thus  demonstrating  that  little  penalty  has  been  paid  in  reducing  the  computa- 
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tional  load.  As  expected,  the  square-root-tree  with  feedback  versions  of  the  algorithm  perform  better 
than  the  basic  version.  There  was  no  discemable  difference  between  the  lattice  version  and  the  array 
version  in  any  of  the  simulations  run  so  far.  This  would  seem  to  demonstrate  the  power  of  the  feedback 
technique  in  improving  the  numerical  accuracy  of  these  algorithms. 


The  relative  effect  of  the  square-root- free  and  the  feedback  techniques  can  be  seen  in  Figure  14. 
This  shows  the  performance  of  the  lattice  algorithms  with  4  bit  mantissas  and  fixed  eigenvalue  spread 
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selling  (W=2.9)  for  the  four  possible  Givens  rotation  algorithms.  From  this  it  can  be  seen  that  there  is 
indeed  a  numerical  advantage  to  avoiding  the  square-root  operation  but  that  the  most  significant  im¬ 
provement  comes  about  by  introducing  the  “error  feedback”. 

All  of  the  above  observations  appear  to  hold  essentially  independently  of  the  eigenvalue  spread. 


8.  CONCLUSION 


We  have  presented  a  derivation,  from  first  principles,  of  a  fast  QRD-based  lattice-ladder  algorithm 
for  solving  the  adaptive  filter  problem  (with  pre-windowed  data).  Contained  within  this  algorithm  is  a 
new  lattice  filter  algorithm  for  least  squares  linear  prediction.  The  extension  to  the  multi-channel  case 
highlights  the  similarity  of  the  adaptive  filtering  algorithm  with  the  more  general  0(p2)  QRD  triangular 
systolic  array.  The  derivation  resented  here  shows  that  other,  recent,  orthogonal  least  squares  lattice  al¬ 
gorithms  are  true  QRD-based  algorithms.  The  relationship  between  these  different  derivations  has  been 
highlighted. 

The  algorithm  has  as  its  root  the  QRD  recursive  least  squares  minimisation  algorithm  and  therefore 
is  expected  to  be  numerically  stable.  Computer  simulations  would  seem  to  confirm  this  expectation:  the 
results  of  the  simulations  of  the  new  algorithm,  using  limited-precision  floating-point  arithmetic,  show 
that  very  little  penalty  has  been  paid  in  reducing  the  computational  load.  The  QRD-based  lattice  algo¬ 
rithm  works  essentially  as  well  as  the  QRD-based  triangular  systolic  array  algorithm  but  requires  only 
O(p)  operations  per  time  instant  as  compared  with  0(p2)  for  the  array. 

Of  the  four  possible  algorithms  for  implementing  the  Givens  rotations,  the  simulation  results  show 
that  the  squarc-root-free  with  feedback  form  of  the  algorithm  is  empirically  better  than  the  standard 
form.  The  former  implementation  coincides  with  the  previously  known  RMGS  lattice  algorithm  of 
Ling.  Manolakis  and  Proakis  [11].  The  algorithm  is  explicitly  presented  in  both  the  normal  (square-root/ 
feedforward)  and  square-root-free  with  feedback  forms. 
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10.  APPENDIX 

10.1.  Joint  Process  Estimation 

It  was  shown  in  sections  4.2  and  4.3  that  a  lattice  filter  algorithm  could  be  developed  to  solve  the 
p-th  order  adaptive  filtering  problem  for  the  previous  time  instant.  Producing  the  adaptive  filtering  re¬ 
sidual  for  time  (n-1)  at  time  n  may  well  be  sufficient  for  some  purposes,  however,  it  is  possible  to  obtain 
the  most  up-to-date  solution  possible  by  including  the  effect  of  the  very  latest  datum  x(n).  Consider  the 
following  data  matrix: 


[xp(“>  yoo  ®„] 


x(l)  £>‘  y(l)  0 

XP-l<n-  Dy(n)  5„_, 


(72) 


where 


Jg'fl)  =  [y(2) . y(n)]‘ 


(73) 


is  the  reference  vector  of  what  may  be  considered  to  be  an  auxiliary  joint  process  estimation  problem. 
Note  that  this  definition  is  more  closely  akin  to  that  for  forward  linear  prediction  (equation  (30))  than 
the  original  definition  of  the  adaptive  filtering  problem  (equation  (6)).  In  equation  (72),  it  has  again  been 
assumed  that  the  data  sequence  x(n)  is  pre-windowed  (i.e.  x(n)  =  0  for  n  <  0). 


By  comparison  with  the  development  of  the  order  update  for  backward  linear  prediction  (section 
4.3)  it  should  be  clear  that  the  matrix  Xp(n),  in  equation  (72),  can  be  triangularised  by  the  series  of  ro¬ 
tations  outlined  in  equations  (47),  (49),  (50)  and  (51).  In  fact 


Qb.p(n)  Qbp(n) 


Qb,p<n-  1)  9 

-  1 

B(n)  ' 

X(l) 

©_ 

o 

_ 1 

p‘  1 

£>  Qp_  ,(n  -  1) 

yf(">  V 

,(n-  l)y(n)  jtn_, 

Rp(n)  up(n)  ap(n) 

°  yp(n)fp(n> 


As  before,  because  certain  matrices  can  be  constructed  directly,  given  the  solutions  to  the  (p-1  )st 
order  problems,  only  the  rotations  summarised  in  equation  (50)  need  actually  be  performed  i.e. 
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10.2  Givens  Rotations 

As  mentioned  in  section  5,  there  are  two  ways  of  calculating  the  parameters  of  a  Givens  rotation: 
with  and  without  a  square-root  operation.  In  this  section  the  basic  mathematics  is  sketched  out. 

A  Givens  rotation  is  a  planar  rotation  that  annihilates  one  component  of  a  2-D  vector.  Assuming 
complex  quantities,  let 


• 

c  s 

Pe(n-l) 

e(nj 

-s  c 

a(n  - 1) 

_  0. 

where,  without  loss  of  generality,  c  is  real  and 

c  fk(n-l)  +  s*  ot(n-l)  =  e(n)  (78) 

c  a(n-l)  -  s  Pe(n-l)  =  0  (79) 

c2+isl2=l  (80) 

now  from  equation  (79)  s  =  c  a(n-l)  /  (k(n-l )  (81) 

then  from  equation  (81),  assuming  £  to  be  real: 

c2(l  +  la(n-l)i2/((3£(n-l))2>=  1  <82i 

hence  c  =  P e(n-l) /  J ($E(n  -  1)) 2  +  |a(n  -  1)  2  (S? i 

s  =  a(n-l)/  n  -  1)) 2  +  |a(n  -  l);2  (84) 

Note  that  e(n)=  V((3c(n-  l))2  +  |a(n-  1)|2  (85) 


In  order  to  avoid  calculating  a  square-root,  we  factorise  the  upper  triangular  matrix  (R  say)  as  fol¬ 
lows: 

R  =  Di/2R  (86) 


40 


where 


0 

0 


(87) 


D  = 


r?,  0 


0  r: 


22  •• 


0  0  ... 


It  is  also  necessary  to  remove  the  factor  of  D1/2  from  any  other  quantities  that  are  operated  upon 
by  the  Q  matrix.  Hence  let 


p£(n  -  1)  P(i(n  -  1)  0 

1 

o 

1 

c 

liw 

1  (I(n  -1)0 

Gf(n)  ap.,(n)  Tp.,(n) 

0  >p-,(n) 

a,(n)  ap_,(n)  1 

so  that 


• 

c  s 

P7e(n  -  1)  0 

-S  c 

L  0  A-.(IV 

1  M(«  -DO 

,(11)  op_,(n)  lj 


jv£(n)  0  |  1  p(n)  v 
!  0  1°  “r(n)  ’j 


(88) 


(89 1 


(where  •  represents  a  quantity  of  no  interest).  Thus  we  find  that 

£(n)  =  P2c(n  -  1 )  +  8p_  ,(n).6tf(n) 2 


(90) 


p(n)  =  (P2c(n  -  l)p(n  -  1 )  +  8p_  ,(n)aj(n)ap_  ^n)) /e(n)  (91) 

=  cp(n -  1)  +  s’ap_  ,(n)  (92) 

where  c=  (P2£(n  -  1)) /£(n)  and  s  =  (5p_  ^nla^n)) /£(n)  <93i 

Similarly  ap(n)  =  ap_  ,(n)  —  a^nljKn  -  1)  (94) 

P2E(n-l)8  .(n) 

using  the  fact  that  Sp(n)  =  - =— - -  =  c8p  ,(n)  (95) 


From  equations  (90)  and  (93)  note  that 
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c  +  Oj(n)s 


1 


(96) 


and  hence,  from  equations  (92)  and  (94),  that 


|i(n)  =  H(n  —  1 )  +  s*  ap(n)  (97) 

This  latter  form  (equation  (97))  of  the  update  for  the  quantity  jl(n)  is  the  most  stable  method  for 
implementing  the  square-root-free  algorithm  and  is  intimately  connected  with  the  "error-feedback"  al¬ 
gorithm  of  Ling  and  Proa Jds[  1 1  ]. 
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10.3.  Algorithms 

These  algorithms,  written  in  pseudo- ALGOL,  calculate  the  adaptive  filter  residual  for  a  p-th  order 
system  fed  with  a  pre-windowed  data  sequence  x(n)  and  a  reference  sequence  y(n).  The  first  algorithm 
uses  the  obvious  implementation  of  Givens  rotations  using  square-roots.  The  second  one  avoids  taking 
square-roots  by  calculating  transformed  quantities  and  implements  the  rotations  via  the  feedback  algo¬ 
rithm. 

10.3.1.  "Normal''  algorithm 

START 

INITIALISE  {all  variables)  :=  0; 

FOR  n  FROM  1  DO 

LET  Of  0(n)  :=  x(n);  o^0(n-l)  :=  x(n-l);  a^n-l)  :=  y(n-l);  y0(n-l)  :=  1; 

FOR  q  FROM  1  TO  p  DO 

LET  ^.,(0-1)  :=  M,,,- 1(n-2))2  +  |abq.1(n-l)|2; 

IF  £b,q.i(n-l )  =  0  THEN  LET  cf  q  :=  1 ;  sf  q  :=  0 

ELSE  LET  cfq  :*  peb  q.,(n-2)  / eb  q.,(n-i);  sfq  :=  04,  q.,(n-l)  / £b,q.,(n-l) 

ENDJF; 

LET  Mf  q-i(n)  :=cf  qppf  q.1(n-l)+  sj  ^ q.,(n); 

«f,q(")  :=  cf,q  °f,q-l(n)  '  sf,q 
Pq_,(n-1)  :=cfqPpq.,(n-2)+  s*  qaq.,(n-l); 
aq(n-l)  :=  cf  q  aq.,(n-l)  -  sf  q  ppq.,(n-2); 
yq(n-l);*c^yq.,(n-l); 

COMMENT  prediction  residual  ef  p(n,n)  =  yq(n- 1 )  04  q(n)  COMMENT 

ep(n-l  ji-1)  =  yq(n-l)  aq(n-l)  COMMENT  q-th  order  filtered  residual  COMMENT 

LET  £f,q-j(n) :«  ,(n-  1))2  + jcxf  q_,(n)p; 

IF  ^.,(0)  =  OTHEN  LET  cbq  :=  1;  sbq  :=  0 

ELSE  LET cb  q  :=P£f,q.i{n-l)/£fq.,(n) ;  sbq  :=af.q-i(n)/Ef.q-l(n) 

END_IF; 

LET  ^b.q-l<n  l)  -cb,qP^b,q-l<n-2>4  sb.qab.q-l<n-,>; 

«bfl(")  '=  cb,q  ab,q-l(n',>  -  sb,q  K.q-l(n‘2); 

COMMENT  Yq(n)  :=  cbq  Yq_i(n-1 );  backward  prediction  residual  eb  p(n.n)  :=  YqCnl  cx^  qfn) 
COMMENT 
END_DO 
END_DO 
FINISH 
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10.3.2.  “Square-root-free  with  feedback"  algorithm 

START 

INITIALISE  (all  variables)  :=0; 

FOR  n  FROM  1  DO 

LET ef  0(nji-l)  :=  x(n);  eb  0(n-l,n-2)  :=  x(n-l);  eo(n-l.n-2)  :=  y(n-l);  80(n-l)  :=  1; 

FOR  q  FROM  1  TO  p  DO 

LET  Eb,q-i<n-l>  :=  P2eb>q.,(n-2)  +  8q.,(n-l)  lebq.,(n-l,n-2)l2; 

IF  eb q.,(n-l)  =  0 THEN  LET  Cf  q  :=  1;  sf  q  :=  0 

ELSE  LET  cf  q  :=  P2£b.q.,(n-2)/  Eb >q.,(n-l);  sfq  :=  5q.1(n-l)eb  q.1(n-l  ji-2)  /  Ebq.,(n-1) 
END_IF; 

LET ef  q(n,n-l)  :=  ef  q.,(n,n-l) -  eb q.,(n-U-2)  |Ifq.,(n-l); 

Pf,q-l<n) :=  ff,q-l(nl) +  s’  qef  q(njl-l); 

eq(n- l^i-2)  :=  eq.,(n-l,n-2)  -  eb  q.,(n-l.n-2)  iIq.,(n-2); 

Pq-l(n-l):=  4q.i(n-l)  +  sj  qeq(n-l,n-2); 

6q(n-l):=Cfq8q.,(n-l); 

COMMENT  prediction  residual  ef  p(n,n)  =  5q(n-l)  ef  qCnji-l)  COMMENT 

eq(n- 1  ji- 1 )  :=  8q(n- 1 )  eq(n- 1  ji-2)  COMMENT  q-th  order  filtered  residual  COMMENT 

LET  £f  q.,(n)  :=  P2ef  g.^n-l)  +  8q.,(n-l)  lefiq.,(n.n-l)l2; 

IF  ef  q.,(n)  =  OTHEN  LET  cbq  :=  1;  sb  q  :=  0 

ELSE  LET  cbq  :=P2E/,q.,(n-l)/£fq.,(n) :  sbq -S^ffn-Defq.jtn.n-D/E^.jln) 
END_IF: 

LET ebq(n,n-l)  :=  ebq.,(n-l.n-2)  -  ef  q.,(n.n-l)  4biq.,(n-2); 

^b,q-l^n"'-* ,=  ^b,q-l^n'^^+  ^b. qeb,q-l^n,n‘^: 

COMMENT  8q(n)  :=  cbq  8q_j(n-l);  backward  prediction  residual  ebp(n,n)  :=  8q(n)ebq(n,n  l ) 
COMMENT 
END_DO 
END_DO 
FINISH 
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10.4  Post-processor  Architectures 


The  algorithm  presented  in  this  memorandum  solves  a  canonical  adaptive  filtering  problem  where 
the  reference  signal  (y(n))  is  known.  There  are  situations,  however,  when  no  such  reference  signal  is 
available  and  a  constrained  least  squares  minimisation  is  required  (e  g.  wide-band  MVDR  beamform¬ 
ing).  The  algorithm  also  calculates  the  adaptive  filtering  residual  directly  without  finding  the  optimum 
coefficients.  It  is  useful  in  system  identification  problems  to  know  what  the  optimum  coefficients  arc 
and  as  such  the  algorithm  derived  so  far  cannot  be  used.  All  of  the  above  draw-backs  also  apply  to  the 
triangular  systolic  array  (28]  and  a  series  of  pre-  and  post -processor  architectures  have  been  developed 
to  over  come  them  [15](25][26](30].  It  is  of  great  interest,  therefore,  to  consider  the  application  of  these 
techniques  to  enhance  the  lattice  algorithm  presented  here. 


10.4.1  Parallel  Weight  Extraction 

It  has  been  shown  (26]  that  it  is  possible  to  extract  the  optimum  weight  vector  (e.g.  to  in  equation 
(4))  in  parallel  with  the  computation  of  the  standard  triangular  systolic  array  (Figure  16).  The  rotation 
parameters  (Qp(n))  are  passed  across  into  an  "inverted”  triangular  array  of  processors  that  apply  these 
rotations.  The  data  that  is  fed  into  this  array  is  a  vector  of  zeros.  As  this  vector  passes  down  the  array  it 
is  transformed  into  the  optimum  weight  vector. 


Now  if  the  data  x(n)  is  in  fact  delayed  samples  from  a  time  series,  then  the  system  is  an  adaptive 
filter  and  the  left-hand  triangle  in  Figure  16  can  be  replaced  by  a  lattice  filter.  Note  however  that  the 
rotation  parameteis  (Qp(n))  are  still  produced  by  the  lattice.  This  can  be  seen  as  follows,  consider  the 
series  of  Givens  rotations  that  constitute  the  matrix  Qp(n).  Recall  that  these  operations  are  intended  to 
annihilate  the  new  data  vector  by  rotation  against  the  previous  version  of  the  triangular  matrix  Rp(n-1 ) 
(see  equations  (22)  and  (23)): 
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The  sequence  of  Givens  rotations  are  constructed  15]  such  that  the  element  x(n)  is  annihilated  first, 
by  rotation  against  the  top  row  of  the  matrix  pR^fn-l),  next  the  second  element  of  the  bottom  row 
(x(n-l))  is  annihilated  by  rotation  against  the  second  row  of  (JRp(n-l)  and  so  on.  Coasidcr  the  situation 
after  the  first  q  <  p  (say)  elements  of  the  bottom  row  have  been  annihilated: 


QqQq-1"  Ql 


Qp(n  -  1)  £> 


B(n)Xp(n) 


Rq(n)  Ub.qfr) 

O  P 

O  P 

.  “b.qt") 


s 

T 

T 

O 


(99) 


where  Qn  is  a  Givens  rotation  and  we  have  used  the  fact  that  the  (q+ 1  )st  diagonal  element  of  Rp(n- 1 )  is 
eb,q(n-l)  (see  last  paragraph  of  section  4.2).  The  identity  of  the  elements  ji^n)  and  Op^fn)  follows  from 
the  underlying  data  matrix  and  the  fact  that  upper-triangular  matrix  R^(n)  has  been  formed.  It  is  clear 
that  the  Givens  rotation  Q,,.,|  is  that  which  annihilates  the  elemem  a^fn)  by  rotation  against  (kbq(n-l ) 
but  this  is  exactly  the  definition  of  Qfs,|(n+1)  (see  equation  (41)).  Thus  the  lattice  algorithm  does  cal¬ 
culate  the  rotations  that  constitute  Qp(n)  and  it  is  possible  to  “bolt-on"  the  weight-extraction  processor 
to  the  lattice  (Figure  17). 
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10.4.2  Linear  constraints 

It  has  been  shown  [25]  how  the  triangular  systolic  array  can  be  adapted  to  solve  a  linearly  con¬ 
strained  least  squares  problem.  The  technique  is  to  use  a  trapezoidal  preprocessor  in  front  of  the  trian¬ 
gular  array  (Figure  18).  Briefly,  system  works  as  follows:  the  constraint  vectors  are  first  fed  into  the 
combined  array  followed  by  the  data  in  the  usual  way  (note  that  the  combined  array  is  triangular).  After 
the  constraint  vectors  have  been  entered  into  tire  array,  the  preprocessor  section  is  "frozen’-,  that  is  to 
say  the  stored  parameters  that  are  normally  updated  from  time  instant  to  time  instant  are  kept  constant. 
As  the  data  subsequently  passes  through  the  preprocessor,  the  constraints  are  incorporated  into  the  data 
before  the  least  squares  minimisation  is  performed  in  the  standard  triangular  array. 

It  would  be  very  useful  to  be  able  to  solve  linearly  constrained  adaptive  filtering  problems  effi¬ 
ciently.  However,  when  a  constrained  least  squares  minimisation  problem  is  re  formulated  as  a  canon¬ 
ical  one  (see  [25])  the  data  matrix  does  not  have  Toeplitz  symmetry  (even  if  the  input  data  does).  Thus 
it  is  not  possible  to  “collapse”  the  triangular  array  into  a  lattice.  It  would  be  interesting  to  know  what 
function  the  lattice  equivalent  to  the  architecture  shown  in  Figure  18  performs.  By  equivalent  architec¬ 
ture  we  mean  a  lattice  where  the  first  q  (say)  lattice  stages  were  operated  frozen  mode  once  the  con¬ 
straint  data  had  been  fed  into  them. 

The  pre-processor  architecture  for  implementing  the  linear  constraints  has  many  desirable  features 
however  it  is  not  the  only  method  for  achieving  this  aim  [30],  It  is  possible  to  apply  constraints  to  the 
least  squares  minimisation  problem  using  a  post-processor:  this  architecture  is  a  modification  of  the 
systolic  QRD-based  MVDR  beamforming  method  [15]. 
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The  MVDR  beamlorming  problem  is  a  least  squares  minimisation  problem  with  a  single  linear 
constraint:  namely  that  the  antenna  gain  in  the  “look"  direction  is  fixed.  The  QRD-based  architecture 
she  xrn  in  Figure  19  can  be  used  to  produce  a  MVDR  residual.  The  triangular  array  is  fed  with  data  to 
initialise  it,  then  the  constraint  vector  is  fed  into  the  array  which,  operating  in  the  frozen  mode,  outputs 
a  vector  that  is  loaded  in  to  the  new  array  on  the  right  hand  side.  The  -riangular  array  is  then  return  to 
its  normal  adaptive  mode  and.  as  more  data  arrives,  it  passes  the  rotation  information  Qp(n)  to  the  post¬ 
processor.  The  output  of  the  combined  structure  is  the  MVDR  beamformer  residual. 

Geariy  because  the  QRD-based  lattice  is  "black  box"  equivalent  to  a  triangular  array  fed  by  a 
tapped  delay  line,  we  may  use  the  same  technique  to  implement  a  MVDR  spectrum  analyser.  The  main 
advantage  to  this  post-processor  MVDR  technique  is  that  several  single  constraints,  or  look  directions, 
can  be  implemented  simultaneously  since  the  rotation  parameters  in  Qp(n)  can  be  fed  to  any  number 
of  post-processors.  The  main  draw-back  is  that  the  magnitude  of  the  numbers  stored  in  the  post-proces¬ 
sor  continually  diminishes  as  lime  progresses  and  the  processor  requires  periodic  re-initialisation  (as 
above)  if  serious  numerical  problems  are  to  be  avoided. 

Yang  and  BOhme  (30)  have  shown  that  by  combining  the  outputs  of  several  different  MVDR  prob¬ 
lems.  the  solution  to  a  multiply  constrained  problem  can  be  generated  (Figure  20).  Apart  from  the  prob¬ 
lem  with  the  MVDR  post-processor  mentioned  above,  this  technique  has  the  draw-back  that  the  final 
trapezoidal  array  has  to  implement  a  hyperbolic  rotation  which  is  numerically  unstable. 
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11.  ANNEX 


As  was  mentioned  in  section  7,  many  computer  simulation  experiments  were  undertaken  and  only 
the  salient  points  have  been  discussed  so  far.  For  the  sake  of  completeness,  all  of  the  simulation  results 
run  to  date  are  included  in  this  annex. 

Each  of  the  following  graphs  has  a  key  which  lists  all  the  relevant  parameters.  The  following  is  a 
short  explanation  of  each  item: 


FILTER  ORDER . the  order  of  the  adaptive  filter  (always  =11). 

CHANNEL  ORDER . the  order  of  the  “raised  cosine”  channel  (always  =  3). 

REF.  DELAY . the  amount  of  delay  in  the  reference  signal  path  (always  =  7). 

APRIORI  ERROR . confirmation  of  the  type  of  error  being  plotted. 

WINDOW . value  of  |3  (see  equation  (2)). 

PRE-LOAD . value  used  to  initialise  energy  terms. 

BITS . number  of  bits  in  the  mantissa  (0  =  double  precision). 

ZERO . value  below  which  quantities  are  considered  to  be  effectively 

zero. 

ITERATIONS . number  of  realisations  of  the  simulation  used  to  calculate  the 

ensemble-average  error. 

SPREAD . value  of  W  (see  equation  (71)). 

VARIANCE . variance  of  gaussian  measurement  noise  (always  =  10‘3). 

G.  APPROX . number  of  iid  random  variables  added  together  to  form  a 

gaussian  random  variable  (central  limit  theorem!). 

SEED . seed  value  for  random  number  generator. 

Several  different  algorithms  were  used  in  the  simulations  and  each  one  was  allocated  a  different  tag: 

NORMAL  LATTICE . Lattice  algorithm,  normal  Givens  rotations. 

SQROOT  LATTICE . Lattice  algorithm,  square-root-free  with  feedback  Givens 

rotations. 

E/FB  LATTICE . Lattice  algorithm,  square-root  with  feedback  Givens 

rotations. 

SF-EF  LATTICE . Lattice  algorithm,  square-root-free  without  feedback  Givens 

rotations. 

FULL  SQRT  ARRAY . Lattice  algorithm  (sic),  normal  Givens  rotations  but  with 

double  precision  arithmetic  inside  the  square-root  operator. 

NORMAL  ARRAY . Triangular  systolic  array  algorithm,  normal  Givens  rotations. 

SQROOT  ARRAY . Triangular  systolic  array  algorithm,  square-root-free  with 

feedback  Givens  rotations. 
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